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SOLVING  ALGEBRAICALLY  EXPLICIT  DAES 
WITH  THE  MANPAK-MANIFOLD-ALGORITHMS* 


BY 

Werner  C.  Rheinboldt^ 

Abstract.  Recently,  the  author  introduced  a  package  of  algorithms,  called  M.\NPAK,  for 
effective  computations  on  implicitly  defined  submanifolds  of  R" .  Here  algebraically  explicit 
differential  algebraic  equations  (DAEs)  are  considered;  that  is,  D.4Es  in  which  either  the 
algebraic  equations  and/or  the  algebraic  variables  are  explicitly  specified.  Existence  proofs 
for  several  types  of  such  DAEs  of  index  one.  two,  and  three  are  given  which  directly  suggest 
computational  approaches.  This  is  used  in  the  development  of  solution  algorithms  for  these 
DAEs,  all  of  them  intrinsically  based  on  the  MANPAK  routines.  Some  numerical  examples 
for  the  methods  are  included. 


1.  Introduction. 

In  [Rh96)  a  package  of  algorithms,  called  MANPAK,  for  computations  on  implicitly 
defined  submanifolds  of  R”  was  presented.  More  specifically,  we  work  with  a  sufficiently 
smooth  mapping  F  :  R”  R^,  d  =  n  —  m  >  0,  which  is  a  submersion  on  its  zero 
set  Af  =  F“^(0),  whence  Af  is  a  d-dimensional  submanifold  of  R”.  MANPAK  includes 
algorithms  for  computing  local  parametrizations  on  such  an  implicitly  defined  manifold  A/ 
and  its  tangent-bundle  TA/,  bs  well  as  other  useful  quantities  on  A/  including  sensitivity 
measures  and  the  second  fundamental  tensor. 

In  [Rh96)  several  applications  of  the  MANPAK  algorithms  were  mentioned.  Here  we 
consider  another  such  application,  namely  the  solution  of  certain  algebraically  explicit 
differential  algebraic  equations  (DAEs);  that  is,  of  DAEs  in  which  either  the  algebraic 
equations  and/or  the  algebraic  variables  are  explicitly  specified.  Algorithms  are  described 
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for  solving  six  types  of  DAEs  of  index  one  to  three,  all  of  them  are  based  intrinsically  on  the 
MANPAK  routines.  These  solvers  have  been  implemented  as  a  collection  of  FORTRAN-  /  / 
subroutines.  As  with  all  the  MANPAK  methods,  the  solvers  are  intended  for  application 
to  small  or  medium-sized  DAE  problems,  mainly  because  they  involve  many  dense  matrix 
computations. 

Section  2  presents  existence  theories  for  algebraically  explicit  DAEs  in  a  form  exhibiting 
the  algorithmic  approach.  Then,  in  Section  3,  the  computational  algorithms  for  these  t^ies 
of  DAEs  are  developed.  Finally,  Section  4  gives  some  numerical  results  obtained  with  the 
implementations. 


2.  Existence  Results. 

We  begin  with  autonomous  DAEs  of  the  form 

Fi(u.  u',iy)  =  0 
F2(u)  =  0, 


iN2) 


under  the  following  condition: 

Assumption  N2;  Let  >  0,  >  0,  be  integers  such  that  ku  -I-  b«r  =  mi  -l-m2, 

ku  >  m2,  and  let  C  C  R''*  be  nonempty,  open  sets.  Assume  that  (i) 

Fi  :  El  R"*'  is  C*  on  Ei  =  Eu  x  R^*  x  OV  E2  :  Eu  >-*  R*"’  is  on  Eu,  (m) 
ranJf  DF2(u,<)  =  m2  on  M2  =  Fj-'CO)  =  {u  €  E„;F2(u)  =  0}  whence  M,  is  a  C*  - 
submanifold  of  R*“  of  dimension  d  =  k„  —  m‘2;  (iv)  the  n  x  n  matrix,  n  =  ku  +  k^, 


(2.1) 


DpF\{u,p,w)  Du)F\{u,p,w) 


DF^iu)  0 

is  nonsingular  on  the  set  K  =  {{u,p,w)  €  E\;  Fi{u,p,w)  —  0,(u,p)  €  rM2}  where  TM2 
is  the  tangent  bundle  of  M2. 

We  are  interested  in  a  solution 


nj’H 


u  ;  Eu.  Eu’^  t  €  cf , 


(2.2) 


of  (N2)  which  is  C  on  an  open  inter\-al  J  containing  the  origin  and  satisfies  the  initial 
condition 

(2.3)  u(0)  =  tin-  ii'iO)  =  Po-  ^'^(0)  =  ’^’0- 

For  some  existence  results  on  such  initial  value  problems  we  refer  to  iIUiOl],  [RaFUiOl],  and 
[RaRh94a).  Here  a  proof  is  given  that  exhibits  directly  a  computational  approach. 

Under  Assumption  N2  the  index  of  (N21  equals  two  if  >  0  and  >  0,  one  if  either 
=  0  or  m2  =  0,  and  zero  if  k,„  =  0  and  =  0.  This  is  reflected  in  the  label  (N2) 
where  “N”  stands  for  “nonlinear"  and  “2"  indicates  the  (generic)  index.  We  shall  assume 
first  that  m2  >  0. 

Given  (uo,po,uto)  €  /i,  let  (V^cp)  be  a  local  parametrization  of  near  uo:  that 
is,  C  R**  is  an  open  neighborhood  of  the  origin,  and  the  C* -mapping 

(2.4)  (p  :  V**  R*''’' ,  ^{y)  €  M2,  Vy  €  V'^,  ^(0)  =  «o 

is  a  homeomorphism  of  V*  onto  its  image  and  an  immersion  on  V*.  Then  ( xR'^, 

is  a  C'  loctd  peirametrization  of  TM2  near  (uo,po)-  In  particular,  we  have  po  =  ^>^(0)-o 

for  some  zo  €  R**- 

With  fixed  y  €  V*  consider  the  nonlinear  system  of  equations 

(2.5)  H{y\ z,  u>)  s  F^ (^(y),  D<^{y)z,  w)  =  0, 

in  the  unknown  z,iu,  for  which  ff(0:  zo.u'o)  =  0.  The  Jacobian  of  H  is 

(2.6)  DH{y\2,w)  =  ( DpFi(<p(y). T>v;(y)z.  u--)I?<Pi(0)  DuS^{^{y)^D^{y)z.xv)) . 

and  i?if(0;  2o,tyo)  is  nonsingultu:.  In  fact,  if  (^1,^2)  €  R**  x  R^""  is  a  null-vector  of  the 
matrix  (2.6)  then  {Dip{y)^i,i2)  is  a  null-vector  of  (2.1)  whence  I>^(y)^i  =  0,  6  =  0- 
auid  therefore  also  =  0  because  of  the  injectivity  of  D^piy).  Thus,  the  implicit  function 
theorem  applies  to  the  system  (2.5).  This  proves  that,  after  shrinking  V'  if  needed,  there 
exists  an  open  neighborhood  Vi  €  R^  x  E^.  and  a  unique  -mapping 

Xj*  ;  if  c  R*^  X  R*^”'.  'l'(y)  =  (t’i(y),tt’2(y))  €  w.  Vy  €  V*'. 

:5 


(2.7) 


such  that,  for  any  y  €  V**,  the  unique  solution  of  (2.5)  in  U  is  gi\'en  by  (r,u>)  —  ^'(y). 
Consider  the  local  initial  value  problem 

{LN2)  y'-Mvh  y(0)  =  o, 

which,  by  standard  ODE-theory,  has  a  unique  C^-solution  y  :  J  V*  on  an  open  interval 
J"  containing  the  origin.  Then,  by  construction  of  0  definition  of  (p  we  have 

which  shows  that  u(t)  =  <^(y(0)»  ^  solution  of  (N2).  Evidently, 

this  solution  satisfies  the  initial  conditions  (2.3).  In  other  words,  we  proved  that,  for 
(tto,po»«>o)  €  K,  there  exists  a  local  C*  solution  of  the  initial  value  problem. 

Conversely,  let  (2.2)  be  any  C* -solution  of  (N2)  on  an  open  interval  J  containing  the 
origin.  Then  we  have  Fi(tt(<),tt'(<),  «;(<))  *  0  and  f2(u(<))  =  0  for  ail  t  €  */•  Hence 
by  differentiation  with  respect  to  f  we  see  that  DF2{u{t))u'it)  ^  0  for  all  t  €  ,/  and 
therefore  (u(t),u'(t),u>(t))  €  /f  for  all  t  €  J.  Thus,  the  above  results  apply  at  any  point 
of  this  solution.  Since  (N2)  is  autonomous,  it  suffices  to  work  with  the  point  (uo,po,  u'o)  = 
(u(0),u'(0),w(0)).  Then,  with  a  local  parametrization  (V*^,!,?)  of  Afj  at  uo  €  A/j  and 
zo  €  such  that  D(p{0)zo  a=po,  the  -  mapping  'I'  of  (2.7)  is  well  defined  on  V**.  For 
the  local  curve  y  ^  V^,  y  **  ou,  we  have  «(<)  =  9(y(0)*  *^^(0  ~  •^‘^(y(O)y  (0» 
y(0)  «  0  which  shows  that  Fi(<p{y{t)),Dtp{yii))y'it),wit))  =  0,  for  aU  <  €  J.  In  view  of 
the  uniqueness  of  the  solution  of  the  equation  (2.3),  this  implies  that  ’i’(y)  =  (y*(t)i  w(<)) 
and,  hence,  that  the  local  curve  is  a  solution  of  the  initial  value  problem  (LN2). 

Altogether  this  proves  the  following  result: 

Theorem  2.1:  Suppose  that  Assumption  N2  bolds  for  the  DAE  (N2).  Then,  for  any  point 
(uo,Po>  Wo)  €  K  there  exists  a  unique  -solution  (2.2)  of  (N2)  on  some  open  inten'al  J. 
to  €  JT  which  satisfies  the  initial  conditions  (2.3). 

As  noted  earlier,  the  index  of  (N2)  reduces  to  one  if  either  Tn2  =  0  or  fcu,  =  0.  Both 
cases  are  covered  by  the  theorem.  But,  for  m2  =  0  the  manifold  M2  equals  the  open  set  Eu 
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and  the  proof  of  Theorem  2.1  no  longer  suggests  a  very  efficient  computational  approach. 
.4ccordingiy,  we  give  here  a  separate  proof  for  this  special  case. 


For  mo  =  0,  we  write  (N2)  in  the  form 


(Nl) 


\  J’(u,p,u;)  =0 


where,  in  analogy  to  Assmption  (N2),  we  now  use  the  condition: 


Assumption  Nl;  Let  /;«,  >  0,  m  =  !;«  +  ky,,  be  integers  and  Ey  C  R*“,  Ey,  C  R 

nonempt}',  open  sets.  Assume  that  (i)  F  :  Ei  R”*  is  on  E  —  Ey  x  Ey  'x  Ey,,  (ii) 
themxm  matrix  ( DpF(u,p,w)  DwF(u,p,  w) )  is  nonsingular  for  any  (u,p,  w)  €  M  = 
F~^(0),  whence  M  is  a  submanifold  of  R”,  n  =  2  *  ky  -h  ky,,  of  dimension  d  =  ky. 


Let  (V^,9)  be  a  C*-local  parametrization  of  M  near  some  point  (uo,poi  u.’o)  €  ilf ,  where 
we  use  for  9  the  component-notation 

9  :  V**  R"  S  R*“  X  R**  X  R^’" 

9(1/)  =  (9i(v),92(y)*<i^’3(y))  €  M,Vy  €  V**,  9(0)  =  (no,po.u'o)» 

For  given  y  €  V*  the  dxd  matrix  I?9i(y)  is  nonsingular.  In  fact,  suppose  that  D9i(y)^  =  0 
for  some  ^  €  R**.  Then 

0  =  DF{^{y))D<piy)f,  *  ( DpF{,ipiy))  Dy>F{<p{y)) ) 

and,  by  condition  (ii)  of  Asstimption  Nl,  the  matrix  on  the  right  side  is  nonsingular  whence 
D<^2{y)i  =  0,  and  D93(y)^  =  0.  Thus  altogether  we  have  D<^{y)i  =  0  which  requires  that 
^  =  0. 

Consider  now  the  local  initial  value  problem 


U^3{y)J 


(liVl)  D<piiy)y'  =  92(y),  y(0)  =  0, 

which  by  ODE  theory  has  a  unique  C* -solution  y  J  V**  on  some  open  interval  J 
containing  the  origin.  Then,  by  construction,  we  have 

0  =  F(9(y))  =  F(9i(y(<)),  {y{t))y\t\ <P3(y(t))h  Vt  €  J 


which  shows  that  u(t)  =  v?i(y(t)),  w{t)  =  <^3(i/(0)>  define  a  C’  solution  of  (Nl).  Moreover, 
because  of  i^(0)  =  (uoiPoi  11^0)1  the  initial  conditions  (2.3)  hold. 

Conversely  let  (2.2)  be  any  C* -solution  of  (Nl)  on  an  open  interval  J  containing  the 
origin.  Then,  clearly,  {u{t),u'{t),wit))  €  M  for  all  t  €  J  and  the  earlier  considerations 
apply  at  any  point  of  this  solution.,  .^s  before,  it  suffices  to  work  only  with  the  point 
(uo,po,u>o)  =  (u(0),tt'(0),u>(0)).  Let  (V‘',<p)  be  a  C*  local  parametrization  of  M  at  this 
point  and  consider  the  locad  curve  y  :  J  ^  y  =  (Pi"’  °  ”  9(v(0)  and 

u'(f)  =  92(y(t))  imply  that  D<p{y(t))y'{t)  =  92(y(0)  for  t  €  J  and  from  (p(0)  *  uo  = 
u(0)  =  <p(y(0))  it  follows  that  y(0)  =  0.  Hence,  y  is  a  solution  of  the  initial  value  problem 

(LNl). 

With  this  we  obtained  the  following  result: 

Theorem  2.2:  Suppose  that  Assumption  Nl  holds  for  the  DAE  (Nl).  Then  for  any  point 
(uQfPo,u)o,to)  €  M  there  exists  a  unique  C' -solution  (2.8)  of  (Nl)  on  some  open  interval 
J,0  €  J,  which  satisfies  the  initial  condition  (2.7). 


These  results  extend  also  to  second  order,  algebraically  explicit  DAEs  which,  in  general, 
have  index  three.  As  illustration,  we  consider  here  only  the  speciztl  case  of  quasilinear 
systems  of  the  form 

f  A(u, u')u"  +  B(u,  u')w  =  G{u,  u') 

(Q3)  i 

I  F(u)  =  0 

under  the  following  condition: 


Assumption  Q3:  Let  ku,mj  >  0,  k^,m2  >  0,  be  integers  such  that  —  mj  +m2, 

ku  >  m2,  and  £„  C  £u»  C  nonempty,  open  sets.  Assume  that  (i)  A  :  Ej 

and  B  :  El  ^  ),  G  :  Ei  ^  R’”'  are  C*  on  £j  =  £,  x  R*«;  (ii) 

F  :  E2  ^  R"*’  is  on  £2  =  (Hi)  rank  DF(u)  =  m2  on  M  =  £"'(0)  whence  M  is 
a  -submanifold  of  R*“  of  dimension  d  k„  —  m2;  (iv)  the  n  x  n  matrix,  n  =  fc*  + 

f  Aiu,p)  B{u.p)\ 

V££(t/)  0  j’ 


G 


is  nonsinguias  for  any  (u.p)  €  TM. 


Given  (uo,po)  €  TM  let  (V^^)  a  C'*-local  parametrization  of  Mj  near  uq,  and 
(V<<  X  R**,  (9,1)9))  the  induced  local  parametrization  of  TM  near  (uo.po)-  As  before, 
there  is  some  zo  €  such  that  po  =  D9(^^*0' 


We  introduce  the  C* -mapping 


(2.10) 


G  :  V*  X  R*'  R"” 

G(y,z)  =  GMy),D<piy)z)  -  A{^(y),D^{y)z)D^<p{y){z,z).  V(y,r)  €  x  R'', 


and  the  mapping 

f  r  :  X  R**  LiR"") 

(r(!,,;)  =  (.4(9(y),Dv’(y)-1D^j(y)  BMy),D^{y)z)) .  v(y,.-)  €  v' x  R'. 

Then,  the  same  proof  as  for  the  matrix  (2.6)  shows  that  the  matrix  r(y.r)  is  nonsingular 
for  any  fixed  (y,z)  €  V**  x  R**  Let  ■Kd  and  denote  the  canonical  projections  of  the  space 
R'^  X  R**"  onto  its  first  and  second  component,  respectively.  We  consider  the  local  initial 

value  problem 


(1(53)  y'  =  ^,  z' =  7rrfr(y,z)~‘(S(y,r),  y(0)  =  0,  z(0)  =  zq. 

Since  the  right  side  is  of  class  C',  ODE  theory  guarantees  that  (LQ3)  has  a  unique  C 
solution  y  :  J  on  some  open  interv*al  containing  the  origin.  Set 

(2.12)  u(0  =  <p(y(0).  u;(t)  = -„,r(y,z)”‘G'(y(#),^(0),  V<  €  J 


Then  F(u(0)  =  0  for  t  €  <7,  and,  by  differentiation  we  find  that  for  all  f  €  J 
(2.13)  u{t)  =  Dip{y{t))z{t),  u"{t)  =  Dip{y{t))z’{t)  +  DV(y(0)(z(t),z(t)). 


The  second  equation  of  (LQ3)  implies  that 


(  A(u{t).u  {t))D^p{y{t))  B{u{t).u'{t))  ) 


u’(oy 


G(u(t).u'(t))z(n  -  A{u{t),u'it))D'^{yit)){:{t).zit)) 


which,  together  with  (2.13),  shows  that  (2.12)  is  a  C*-soiution  of  (Q3)  for  which 
(2.14)  u{to)  -  uo,  u'ito)  -  po. 

Note  that  the  second  parts  of  equations  (2.12)  and  (2.13)  enforce  initial  \-alues  for  u)(0) 
and  u"(0)  which  therefore  cannot  be  prescribed. 

Conversely,  let  (2.2)  be  any  C^-solution  of  (Q3)  on  an  open  interval  J  containing  the 
origin.  Then  F(u(t))  =  0  and  DFiuit))u'it)  =  0;  that  is,  (u(t),u'(<))  €  TM  for  alU  €  J. 
Thus,  the  previous  consideration  apply  at  any  point  of  this  solution  tind  it  suffices  again  to 
work  with  the  point  (uo,Po,u'o)  =  (u(0),u'(0),u;(0))  where po  =  D<f{0)zo  ior  some  zq  € 
Then  with  a  local  parametrization  ( V**,  ;p)  of  M  at  uo  €  A/,  the  mappings  (2.10)  and  (2.11) 
ane  well-defined  and  the  matrix  r(y,  z)  is  nonsingular  for  smy  fixed  (y,  c)  €  x  For 
the  local  curve  y  '  ^  •— y  =  9~*  ow.  it  follows  that 

*p(y(0)  =*  =  ^^(0*  ■i^<p(y(0)y^^(*)  +  ■i^^<p(y(0)(y*(0»y*(0)  = 

and  y(0)  =  0  and  y'(0)  =  20  •  Moreover,  by  substituting  these  representations  for  u  and  its 
derivatives  into  the  first  equation  of  (Q3)  and  using  (2.10)  and  (2.11).  we  find  that 

/y"(t)\ 

r(y(t),y'(t))  ,  J=G(y(0,ym  vteJ, 

V  w(t)  ) 

whence  y  zuid  2  =  y'  is  a  solution  of  (LQ3). 

This  proves  the  following  result: 

Theorem  2.3:  Suppose  that  Assumption  Q3  holds  for  the  DAE  (Q3).  Then  for  any  point 
(uo,po)  €  El  such  that  (uo,po)  €  TM  there  exist  a  unique  C^-solution  (2.2)  of  (Q3)  on 
some  open  interval  J  C.  Et,  tQ  ^  J  which  satisfies  the  initial  condition  (2.14) 


3.  Computational  Algorithms. 

All  algorithms  in  this  section  (except  one)  are  written  for  non-autonomous  versions  of 
the  DAEs  in  Section  2.  Our  existence  results  are  easily  extended  to  these  non-autonomous 
cases  by  adding,  as  usual,  the  equation  t'  =  0.  The  details  should  not  be  required  but 
some  comments  will  be  provided  with  the  algorithms. 

The  existence  proofs  suggest  that,  in  each  case,  we  should  solve  the  corresponding  local 
ODE  to  compute  the  solution  of  the  (global)  DAE.  Thus  the  process  always  begins  with 
the  construction  of  some  local  parametrization  of  the  manifold  constraining  the  problem, 
followed  by  the  application  of  a  standard  ODE  solver  to  the  resulting  local  system.  When 
the  computed  points  appear  to  leave  the  domain  of  validity  of  the  local  coordinate  system, 
a  new  local  parametrization  is  generated  and  the  process  is  continued  by  applying  the 
ODE  solver  to  the  new  local  ODE. 

For  the  algorithms  discussed  here,  a  new  reverse-communication  version  of  the  explicit 
Runge-Kutta  method  DOPRI5  of  {HNWS7)  was  chosen  as  the  ODE  solver.  Other  inte¬ 
grators  are  readily  applicable  in  this  setting.  For  simplicity,  the  codes  establish  a  neu 
local  parametrization  at  each  accepted  point;  that  is,  after  each  successful  RK-step.  Once 
all  stages  and  the  next  approximate  solution  of  the  local  ODE  have  been  computed,  the 
Rl^-routine  performs  the  standard  error  computation.  If  the  RK-step  is  accepted  then 
from  the  computed  results  a  new  approximate  solution  of  the  DAE  is  determined  and  the 
process  is  repeated  with  this  new  global  point  until  the  required  termination  time  has  been 
reached.  If  the  RK-step  was  rejected,  then,  as  usual,  the  stepsize  is  reduced  and,  from 
the  previous  approximate  global  solution,  the  step  is  repeated  with  the  estimated  smaller 
stepsize.  A  stepsize  reduction  may  also  have  to  be  enforced  if,  during  an  RI\-step  the  local 
parametrization  evaluation  fails  to  converge  which  indicates  that  the  local  vector,  at  that 
stage,  fails  to  belong  to  the  domain  of  the  locad  paLrametrization. 

We  note  that  there  is  no  fimdamental  obstacle  to  retauning  the  local  parametrization 
for  several  successful  RK-steps. 

There  should  be  no  need  to  detail  the  overall  framework  of  the  codes  involving  the  use 
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of  the  ODE  solver.  Instead  we  focus  on  the  principal  aspect  of  the  different  algorithms, 
namely  the  evaluation  of  the  right  side  of  the  local  DAEs.  All  these  algorithms  are  called 
in  two  distinct  modes  to  distinguish  whether  or  not  the  local  parametrization  has  to  be 
established  at  the  particular  point. 

In  the  case  of  (N2)  the  algorithm  is  applied  to  the  DAE 

r  fi(u,u',u;,t)  =  0, 

(4.1)  { 

I  F2(u,t)  =  0 

where  the  equation  =  0  is  not  explicitly  stated.  Note  that  now  A/j  =  ^2  (0)  Is  a 
submanifold  of  R",  (n  =  2k^  +  &„,  + 1),  of  dimension  d  =  /:«  + 1.  Accordingly,  for  the  local 
paremietrization  (^*,1,5)  at  smy  (uoito)  €  Afj  the  component  notation 

(4.2)  9  :  V*  R*“  X  R\  ^(y)  s  (v>i(y),V’j(y))  €  x  .S(,  Vy  €  V**,  (uo,to)  =  ^>(0) 


is  introduced.  The  evaluation  of  the  right  side  of  the  local  ODE  then  requires  the  solution 
of  the  nonlinear  system 


(4.3) 


Hiy,z,w)  = 


__  f  Fi(9i(y),D(,?i(y)r,u?,<r’2(y))\  _ 


I 


Ds9j(y)r  -  1 


=  0. 


With  the  solution  function  of  (2.7)  we  have  here  Dv?2(y)ti'i(y)  =  1  for  all  y  €  V**  which 
implies  that  now  the  relation  between  the  local  and  global  solution  is  u(0  =  <r’(y(t  “  to)), 
w(t)  =  V'2(y(<  -  to)). 


For  the  solution  of  (4.3)  a  standard  chord-Newton  method  is  applied  with  the  Jacobian 
of  Af  as  iteration  matrix  typicsJly  evaluated  at  the  time  of  the  establishment  of  the  local 
parametrization;  that  is,  with  a  matrix  of  the  form 


(4.4) 


/{D,Friu,p,w,t)D^M  0) 

V  1)92(0) 


Du,Fiiu,p,w,t)\ 

0  / 


Then  the  algorithm  for  the  evaluation  of  the  right  side  of  the  local  ODE  has  the  form: 


DYN2  Input:  mode 

If  mode  =  ‘init’  then 


10 


Input:  Global  point  (u,p,uj,t); 

Evaluate  DFiu,t)\ 

With  GNBAS  compute  the  basis  of  the  local  parametrization  9  at  (u.t): 

Set  the  local  point  y  =  (0,0)  €  R*“  x 

Set  up  and  factor  the  iteration  matrix  (4.4)  at  the  point  (it,p, uj,  f). 

Use  DGPHI  to  evaluate  D<fi{y)  -  (I>tpi(y),I><,?2(y)) 

Set  the  start  vector  (C»u>)  with  C  = 

Set  the  parametrization  time  tc  =  t; 

Else 

Input:  Local  point  y  in  the  parametrization  (p  at  time  fc,  last  vector  (C/i^^)t 
With  GPHI  evaluate  the  global  point  (u,t)  =  <piy) 

If  GPHI  fails  to  converge  then  force  an  RK-step  reduction 
Use  DGPHI  to  evaluate  D(^{y)  —  {D(p\{xj),D‘P2{y)) 

Set  the  start  vector  (C,u;)  —  (OttVf, 

End  If; 

Use  the  chord- Newton  method  with  the  current  iteration  matrix  to  solve 
the  nonlinear  system  ff(0,  C,u>)  =  0  of  (2.5)  for 
If  the  method  diverges  Then  error  return 
Set  p  =  D<pi(y)(C); 

Output:  y'  :=  (C,l),  (u,p,t),  ;=  (C.t/^)- 

Since  M  is  here  a  submanifold  of  R*"  x  R^ ,  where  the  one-dimensional  space  represents 
the  t  variable,  the  local  parametrization  is  constructed  by  means  of  the  algorithm  GNBAS 
of  [Rh96]  which  preserves  the  t  variable.  Correspondingly,  for  the  computation  of  the 
derivative  D^{0)  of  the  local  parametrization  the  algorithm  DGPHI  of  [Rh96]  is  applied. 
If,  in  the  ’init’-mode,  we  were  assured  that  (ucPc^^c^ic)  €  K  then  we  would  expect  that 
H{0,z,wc)  =  0  for  r  =  [Dip(0)'^ D<p{0)]~^ D<p{0)'^Pc-  For  the  algorithm  it  was  found  to  be 
more  advantageous  to  enforce  the  validity  of  the  equation  explicitly  by  applying  the  chord 
Newton  method  to  (4.4)  starting  with  given  by  (  =  Dip{0)^{p,  1)  and  w  =  This 
usually  converges  in  a  step  or  two  and  is  no  more  costly  than  the  indicated  direct  evaluation 
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of  r.  Moreover,  no  convergence  indicates  here  that  the  point  is  too  far  away  from  K  to 
be  corrected  onto  that  manifold.  The  algorithm  for  the  solution  of  initial  value  problems 
for  the  nonautonomous  version  (4.1)  of  (N2)  has  been  implemented  as  a  subroutine  suite 
caUed  DAEN2. 


(4.5) 


The  process  simplifies  considerably  in  several  special  cases  of  (N2).  In  particular,  in  the 
quasilinear  case 

A{u,t)u'  +  Biu,t)w  =  Giu,t) 

F2{u.t)  =  0, 

the  chord  Newton  iteration  is  no  longer  any  needed.  Moreover,  only  the  initial  condition 
(uo,to)  €  M2  has  to  given  is  and  no  initied  vectors  po  &iid  wo  are  needed.  For  each  local 
vector  y  €  arising  during  the  Rl^’Step,  we  have  to  set  up  and  solve  the  linear  system 

((AMy))DvM  0)  / 


(4.6) 


I,  Z?v’2(y)  0  / 


=  GMy)) 


to  obtain  z  =  ^i(y)  and  u>  =  r/>2(y).  All  other  parts  of  the  algorithm  remain  the  same. 
The  algorithm  for  solving  initial  value  problems  for  the  special  case  (4.5)  of  (N2)  has  been 
implemented  as  a  subroutine  suite  called  DAEQ2. 

In  the  special  index-one  case  (Nl)  of  (N2)  our  proof  already  indicated  the  considerable 
simplifications  in  the  algorithm.  The  algorithm  is  applied  to 

u'  =p, 

F(u,p,w,t)  =  0, 

where,  again,  the  equation  =  0  is  subsumed.  In  this  case.  Mi  is  a  d  =  -H  l-dinaensional 
submanifold  of  R”,  (n  =  2fc„  +  +  1),  and,  analogous  to  (4.2),  we  write  the  local 

parsunetrization  in  the  component  form 

(4.S)  9  :  V*  «  R",  9(y)  =  (9l(y),9j(y).^’j(y).>'-4(»))  «  xR*"  xR',  Vy  €  V. 


(4.7) 


Then  the  local  ODE  becomes 
(4.9) 


D‘pi{y)y'  =  <p2{y) 
D(pA(y)y'  =  1 
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and  the  evaluation  of  its  right  side  has  the  following  form: 


DYNl:  Input:  mode 
If  mode  =  ‘init’  then 

Input:  Global  point  (u,p,/); 

Evaluate  DF(u,p,  t); 

With  GNBAS  compute  the  beisis  of  sp  at  (u,p,t); 

Set  the  local  point  y  =  (0, 0)  €  R*“  x  R’ ; 

Parametrization  time  te  —  t; 

Else 

Input:  Local  point  y  in  the  parametrization  (p  at  time  tc 

With  GPHI  evaluate  the  global  point  (u,p,t)  =  (<Pi(y))<r’2(y),‘r?3(y)) 

If  GPHI  fails  to  converge  then  force  an  RK-step  reduction 
End  If; 

With  DGPHI  evaluate  D<p{y)  =  {Dy>i{y),  Dfp2{y)> 

Solve  the  linear  system  D(pi(y)C  =  (Pily)  for  Ci 
Output:  y'  :=  (C.l),  (u,P,0- 

The  algorithm  for  the  solution  of  initial  value  problems  for  the  nonautonomous  version 
(4.7)  of  (Nl)  has  been  implemented  as  subroutine  suite  called  DAENl. 

In  analogy  to  (Q2)  we  may  admit  also  a  specialization  of  (Nl)  to  quasilinear  form.  Such 
index-one  systems,  albeit  in  the  autonomous  form 

(  A{u)u'  =  G{u) 

(4.10)  { 

I  F2(u)  =  0, 

have  been  considered  earlier  in  [RaRh94c]  in  connection  with  the  computation  of  impasse 
points  (see  also  (RaRh94b]).  In  this  case  the  reduced  system  has  the  form 

(4.11)  .4((p(j/))Fip(y)y'  =  G(ip(y)) 

and  it  is  obvious  how  to  change  the  above  algorithm  DYNl  for  this  case.  The  autonomous 
form  was  retained  here  to  allow  for  a  direct  integration  of  the  solution  algorithm  with 
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the  methods  of  {RaRh94c].  In  our  setting  the  impasse  points  are  those  points  where  the 
reduced  ODE  has  a  standard  singularity,  in  the  sense  that  A{y)  =  A(<p(y))D<,^(y}  satisfies 

dim  ker  A(y)  =  1,  DA(y)(z,  r)  ^  rge  A(y),  V-  €  ker  A(y). 

The  combined  algorithm  has  been  implemented  in  the  form  of  a  suite  of  subroutines  called 
DAESQl  where  “S”  stands  for  “singvdar”  point. 

In  the  case  of  (Q3)  the  algorithm  is  applied  to  the  nonautonmous  version 

f  A(u,u',t)u"  +  B(u,u\i)w  =  G(u,u\t) 

(4.12)  \ 

1  F(u,f)  =  0 

Here  A/j  is  a  submanifold  of  R”,  n  =  2k^  +  fcu,  + 1,  of  dimension  d  =  A 1  and  we  use  the 
component  notation  (4.2)  for  the  local  parametrizations.  The  algorithm  for  the  evaluation 
of  the  right  side  of  the  local  ODE  now  hats  the  form: 


DYQ3  Input:  mode 

If  mode  =  ‘init’  then 

Input:  Globed  point  (u,p,t); 

Evaluate  DF{u,  t); 

With  GNBAS  compute  basis  matrix  Uc  of  (p  at  (u,t); 

Set  the  local  point  y  =  (0,0);  r  =  Ujpe 
Set  the  paramettization  time  te  =  t; 

Else 

Input:  Local  point  y,z  in  the  local  parametrization  <p  at  time  tg 
With  GPHI  evaluate  the  global  point  (u,<)  =  (v’i(y)»  V’2(y)) 

If  GPHI  fails  to  converge  then  force  an  RIC-step  reduction 
End  If; 

Use  DGPHI  to  evaluate  the  derivative  D(p{y)  =  (Dtpiiy),  Dip2(y))', 
Set  p  =  Dipi{y)z; 

Evaluate  q  =  D^<p(y){z,z)\ 

Evaluate  G  —  G(u,p,t)  —  A(u,p,t)q', 
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Evaluate  F  =  (A(u,p,t)Dip{y)B(u,p,t))\ 

Solve  the  linear  system  r(C,uj)^  =  G  for  C  and  iv\ 

Output:  y'  :=  r;  z'  :=  Ct  (u,p,u?,f). 

The  algorithm  for  the  solution  of  initial  value  problems  for  the  nonautonomous  version 
(4.6)  of  (Q3)  has  been  implemented  as  subroutine  suite  called  DAEQ3. 

Evidently,  the  Euler- Lagrange  equations 

'  M(u,t)u”  -t-  DuF{u,t)'^w  =  G{u,u\t) 

(4.13)  < 

F(u,t)  =  0, 

arising  as  models  of  multibody  systems,  are  a  special  case  of  (4.12).  In  that  case,  the 
condition  (iv)  of  Assumption  (Q3)  is  equivalent  with  the  usual  assumption  that  the  mass 
matrix  M  is  definite  on  Tu.tAfj.  In  (RaRh95)  a  different  algorithm  for  the  numerical 
solution  of  (4.13)  was  given  which  uses  the  second  fundamental  tensor  of  the  manifold  M2- 
This  algorithm  has  also  been  added  to  the  present  collection  of  DAE  solvers  under  the 
name  DAEUL3. 

4.  Computational  Examples. 

As  described  in  the  previous  sections,  the  package  of  solvers  for  algebraically  explicit 
DAEs  consists  of  the  following  suites  of  routines: 

(1)  DAENl  for  nonlinear,  index-1  DAEs  (4.1)  of  order  one 

(2)  DAESQl  for  quasilinear,  index- 1  DAEs  (4.10)  of  order  one  with  singular  points 

(3)  DAEN2  for  nonlinear,  index-2  DAEs  (4.2)  of  order  one 

(4)  DAEQ2  for  quasilinear,  index-2  DAEs  (4.5)  of  order  one 

(5)  DAEQ3  for  quasilinear,  index-3  DAEs  (4.12)  of  order  two 

(6)  DAEUL3  for  Euler-Lagrange  equations  (4.13). 

FORTRAN  77  versions  of  the  codes  are  available. 

We  give  here  some  numerical  examples. 
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The  following  index-two  problem  was  given  in  (MT94j 


1 


uj  +  sin(arcos  ui ) - j  u'*  -f- 1 


0 


(5.1) 


Uj  -H  ti?  =  0 

U2  —  logui  =  0 


For  u(0)  =  (1,0),  u;(0)  =  0,  the  exact  solution  is  ul(<)  =  cost,  uJ  =  log  cost,  uj*(t)  =  tant. 
In  [MT94]  the  system  was  integrated  from  t  =  0.5  to  t  =  1.5  using  accurate  starting  values 
2uid,  among  others,  a  BDF-solver  of  order  2  with  fixed  stepsize  10"®  and  a  tolerance  of 
10"®  for  the  Newton  iteration.  The  resulting  global  errors  at  t  =  1.5  were 

|u*  -  ujl  =  0.560(~9),  |U2  -  ujl  =  0.791(-S),  |u;  -  u)*|  =  0.112(-6). 

When  DAEN2  was  applied  to  (6.1)  with  accurate  starting  values  and  a  relative  tolerance 
RTOL  =  10"*,  the  terminal  point  at  <  =  1.5  was  reached  in  38  Rl'C-steps  (no  rejections) 
and  the  final  absolute,  global  errors  were 

\ui  -  ujl  =  2.734(-10),  |U2  “  u;|  =  3.115(-9).|u;  -  u»’I  =  5.476(-8). 

More  illuminating  axe  the  relative,  global  errors 

|u,  -uj|/|uj|  =  3.865(-9),|u2  -uJI/luJI  =  1.176(-9),  |u;  -  u>*|/Kl  =  3.884(-9), 

which  show  no  deterioration  in  the  error  of  the  edgebraic  variable. 

The  second  example  is  the  trajectory-prescribed-path  control  problem  of  [B83]  which 
is  discussed  in  detail  in  [BCP89].  We  refer  to  the  latter  reference  for  the  equations  of  this 
index  two  problem.  It  involves  six  differentizd  variables,  ui, . . .  , ue,  two  algebraic  variables 
1^1,  ^2,  both  of  which  occur  nonlinearly.  The  system  consists  of  six  differential  and  two 
algebraic  equations.  At  t  =  0  the  initial  conditions 

u(0)  =  (0.0, 100000.0, 0.0, 12000.0, -1.0, 7r/4) 

were  imposed  where  U],i(3,U5,us  are  in  radians  Md  U2,U4  in  feet.  The  remaining  initial 


values 


u)(0)  =  (0.046650383, -0.91122917( -3)) 

u'(0)  =  (0.40394369(-3),  -209.42888, 0.40394369( -3),  -34.978260.0.0,0.0) 


where  computed  by  DAENl  and  are  in  agreement  with  those  given  in  (BCP89].  With  a 
relative  tolerance  RTOL  —  10”*  the  system  was  integrated  to  t  =  300.  The  nnai  point 
was  reached  in  44  RI'L-steps,  including  3  rejections,  and,  as  Table  5.1  shows,  agrees  ^er^ 
closely  with  the  reference  values  cited  in  (BCP89). 


DAEN2 

BCP89 

Ui 

0.0727991598693 

0.0727991600 

U2 

14200.8114923 

14200.8114 

U3 

0.0406923052169 

0.0406923053 

U4 

1433.29213943 

1433.29213 

U5 

-0.174532925199 

-0.174532925 

Ufl 

2.35619449019 

2.35619449 

U?i 

0.124888510970 

0.124888511 

U'l 

0.460375007767 

0.460375012 

Table  5.1 

As  an  example  for  DAEQ2  we  use  the  index-two  DAE  proposed  in  [AP91) 

u'l  +  a(t  -  2)u)  *  (a  + 

<x  —  1  . 

(5.2)  -  (a  -  l)ti;  * 

0  =  (2  4-  <)ui  -b  -  4)u2  -  (<^  -b  t  -  2)c‘. 

For  a  =  50  and  the  initial  conditions  u(0)  =  (1,1),  the  exact  solution  is  ui(<)  =  U2(<)  = 
expt,  w  =  (expt)/(t  -  2).  With  a  relative  tolerance  RTOL  =  10"®,  DAEQ2  reached 
t  =  1.0  in  202  RK-steps  (no  rejections)  and  the  maximum  norm  of  the  final,  absolute 
error,  was  0.4(-9).  When  run  with  DAEN2,  (5.2)  happens  to  be  one  of  the  few  problems 
where  it  does  make  a  difference  what  updating  strategy  is  used  for  the  iteration  matrix 
in  the  chord-Newton  process.  When  the  iteration  matrix  is  retained  throughout  an  entire 
RK-step,  then  DAEN2  required  328  RK-steps  (including  97  rejections)  to  reach  t  =  1.0.  On 
the  other  hand,  as  expected,  the  performance  of  DAEN2  was  identical  to  that  of  DAEQ2 
when  the  iteration  matrix  is  updated  for  each  chord-Newton  iteration. 

As  a  first  example  for  DAENl  we  used  the  simple  problem 

(5.2)  u^-b«'*~l=0,  2uu'-w  =  0. 
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For  the  initial  conditions  u(0)  =  0,  u'(0)  =  1.  u;(0)  =  0  it  has  the  exact  solution  u(t)  =  sin  t. 
w{t)  =  sin2t  which  represents  a  horizontal  figure  eight  in  the  plane  and  hence  is  not  a 
submanifold  of  R^.  Note  that  here  the  matrix  (DpFiu.w)  Du>F{u,w))  of  Assumption 
2.2  is  nonsingular  on  My  except  at  the  point  (u,p,w)  —  (1,0,0)  which  turns  out  to  be 
a  removable  singulsirity  of  M\.  In  all  runs,  DAENl  stepped  over  this  point  when  it  was 
reached  for  t  =  (2jt  +  l)7r/2,  it  =  0, 1, . . . .  In  fact,  the  absolute  error,  in  the  maximum 
norm,  remained  about  of  the  same  order  of  magnitude  as  the  chosen  relative  tolerance 
RTOL  =  10“®.  On  the  other  hand,  DASSL  (see  (BCP89))  always  terminated  just  before 
t  =  Tr/2.  A  second  example  is  the  batch  reactor  model  given  in  [BD86] 

u[  +  k2U2W2  =  0 
U2  +  fcltljUfl  --  k-\Wi  +  k2U2W2  —  0 
U3  —  k2U2W2  —  kiUiUt  +  fc-3U>3  =  0 
U4  +  klUiUt  —  k-%W3  =  0 

-  kiU2Ui  +  fc-lU’4  =  0 
Ug  +  kiU2U^  +  k^UiU^  —  k^\W4  —  fc_3U>3  =  0 
Ufl  —  u;i  +  ti;2  +  ti’s  +  ~  a  =  0 

U»2  -  (A'2Ui  )/(/V2  +  U)1 )  =  0 
u>3  —  (/vaUsj/C/Ca  +  tui)  =  0 
tl>4  -  (A'iU5)/(A'i  +  tDi)  =  0, 


where 

it,  =  21.893,  ifc-i  =  2.14(+9),1;2  =  32.318,  Ars  =  21.893,  fc-3  =  1.07(+9), 
Ki  =  7.6o(-18),A'2  =  4.03(-ll),  A'3  =  5.32(-18),a  =  0.0131 


The  initial  conditions  are 

u(0)  =  (1.5776, 8.32, 0.0, 0.0, 0.0, a),  u;(0)  =  (0.79735161(-5),0.79735161(-5),0.0,0.0) 

and  the  values  of  the  derivatives  were  determined  by  the  differential  equations.  For  this 
problem,  DASSL  did  not  start  while  DAENl  reached  t  =  1.0  [hr]  with  229  RK-steps 
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(including  30  rejections).  Problems  of  this  type  suggest  the  need  for  introducing  also  an 
implicit  RK-solver  into  these  DAE-codes. 

Examples  for  use  of  DAESQl  were  given  in  lRaRh94c]  and  for  the  application  of 
DAEUL3  we  refer  to  (RaRh95]  and  the  comparative  study  (RhS95).  DAEQ3  also  runs 
for  these  Euler-Lagrange  examples.  Instead,  we  consider  here  only  the  following  example 
for  DAEQ3  which  is  not  of  Euler-Lagrange  type 


(5.2) 


ui”  4-  uiu?  =  (1  -|-sint)e* 

2  1 

uj”  +  wjty  = 


(1 


0  = 


1  +  t 
2e‘ 


sint 


-h 


eT+c 


(l+<)3  (H.t)2  (1  +  t) 


-  txjUj  +  3uitX2  -  U2C 


For  tho  initial  conditions  n(0)  =  (1,1),  u'(0)  =  (1,-1)  the  problem  has  the  solution 
uj  =  expt,  uj  =  1/(1  + 1),  UJ*  =  sint.  But  in  this  case  the  matrix  (2.11)  of  Assumption 


2.3  has  the  form 


/I  0  «i  \ 

0  1  “2 

\-U?+3u^  -3uiU2  +  6UjU2  -  ujc"*  -  0' 

and,  on  the  exact  solution,  this  matrix  becomes  singular  at  t*  «  0.03756275.  When  started 
at  t  =  0.0,  DAEQ3  produces  a  solution  for  which,  near  t*,  the  algebraic  variable  tends 
to  -oo  while  the  other  variables  remain  close  to  their  exact  values.  As  an  illustration. 
Table  5.2  shows  both  the  solution  obtained  by  the  code  at  t  =  0.037562749184  and  the 


corresponding  exact  solution. 


DAEQ3 

Exact  Sol. 

1.0382771509 

1.0382771461 

U2 

0.96379713463 

0.96379713014 

XV 

-25014116.226 

0.037553916550 

Table  5.2 


The  same  behavior  is  seen  when  the  code  is  run  backwards  from  t  -  0.05  (with  exact 
starting  \'alues),  except  that  now,  as  expected.  w{t)  tends  to  +cx3.  On  the  other  hand. 


19 


when  run  forward  from  t  =  0.05  the  computed  solution  approximates  the  exact  solution 
with  an  absolute  error  of  the  same  magnitude  as  the  given  tolerance. 

Singularities  of  quasilinear.  second  order  ODEs  have  been  studied  in  (L95].  A  discussion 
of  the  relation  between  these  results  and  the  observed  singularity  for  (5.3)  cannot  be  given 
here.  It  is  noteworthy  that  the  singularity  is  not  apparent  in  the  given  DAE. 
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